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ABSTRACT 

We investigate non-linear scaling relations for two-dimensional gravitational collapse in an 
expanding background using a 2D TreePM code and study the strongly non-linear regime 
(£ ^ 200) for power law models. Evolution of these models is found to be scale-invariant in 
all our simulations. We find that the stable clustering limit is not reached, but there is a model 
independent non-linear scaling relation in the asymptotic regime. This confirms results from 
an earlier study which only probed the mildly non-linear regime (£ ^ 40). The correlation 
function in the extremely non-linear regime is a less steep function of scale than reported in 
earlier studies. We show that this is due to coherent transverse motions in massive haloes. We 
also study density profiles and find that the scatter in the inner and outer slopes is large and 
that there is no single universal profile that fits all cases. We find that the difference in typical 
density profiles for different models is smaller than expected from similarity solutions for halo 
profiles and transverse motions induced by substructure are a likely reason for this difference 
being small. 
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1 INTRODUCTION 



Large-scale structures such as clusters of galaxies are very over- 
dense compared to the average density of the universe. It is be- 
lieved that these were formed by the growth of small perturbations 
via gravitational instability (Bernardeau et al., 2002; Padmanabhan, 
2002b; Peacock, 1999; Peebles, 1980). In general, the equations 
that describe the growth of density perturbations due to gravita- 
tional clustering (Peebles, 1974; Peebles, 1980) cannot be solved 
analytically when the over-densities are large and one has to rely 
on N-Body simulations for detailed predictions. Limited informa- 
tion about non-linear gravitational clustering can be obtained by 
using approximations or ansatze, e.g., stable clustering (Peebles, 
1980) for which it is assumed that the infall velocity and Hub- 
ble velocity compensate for each other leading to a stable config- 
uration. Stable clustering allows us to relate the initial spectrum 
of fluctuations with the final asymptotic spectrum of fluctuations 
(Davis & Peebles, 1977). Non-linear scaling relations are a more 
detailed prescription for relating the initial and final spectrum of 
fluctuations (Hamilton et al., 1991; Nityananda & Padmanabhan, 
1994; Peacock & Dodds, 1994; Jain, Mo, & White, 1995; Peacock 
& Dodds, 1996; Padmanabhan, 1996; Padmanabhan et al., 1996; 
Bagla & Padmanabhan, 1997; Bagla, Engineer, & Padmanabhan, 
1998; Smith et al., 2003). The motivation for studying non-linear 
scaling relations is to understand key phases in gravitational clus- 



tering, and to identify the relevant process in each phase (Padman- 
abhan, 1996). 

Existence of non-linear scaling relations implies that gravita- 
tional clustering does not erase memory of initial conditions. This 
provides another useful and a more prosaic motivation for study- 
ing scaling relations, namely to be able to make predictions about 
clustering of galaxies for a given power spectrum (Hamilton et al., 
1991). Scaling relations can also be used to make predictions for 
weak lensing observations in a given model (Kaiser, 1992; Jain 
& Seljak, 1997). These scaling relations are valid for hierarchical 
models where the initial conditions contain fluctuations at all scales 
and the amplitude of fluctuations increases monotonically as we go 
from large-scales to small-scales. 

Non-linear scaling relations indicate that there are three 
prominent regimes in evolution of gravitational clustering. When 
the amplitude of perturbations is small so that the density contrast 
(5 — p/p — 1) is close to zero (|<5| <C 1), mode coupling is not 
important and the evolution closely follows the predictions of lin- 
ear perturbation theory (Peebles, 1974). The power spectrum and 
correlation function evolve without a change in shape in the lin- 
ear regime. As the density contrast grows and becomes comparable 
to unity, motions induced by gravitational collapse start to domi- 
nate over the expansion of the universe. The quasi-linear regime is 
dominated by infall onto density peaks and the correlation function 
grows rapidly in this phase (Padmanabhan, 1996). Gravitational 
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collapse leads to formation of structures in or close to dynamical 
equilibrium and further evolution of density contrast is dominated 
by depletion of average density due to expansion of the universe 
as the density of collapsed structures remains almost constant. We 
shall call this the asymptotic regime. 

Initial conditions with power law power spectra are very use- 
ful for such studies. These power spectra have the form P(k) oc 
k n with n being the index of the power spectrum. We require 
(n + D) > for hierarchical clustering and N-Body simulations 
can be used for only this type of models, indeed even for these 
models there are restrictions on the size of the simulation box in 
order for the simulation to reproduce the model correctly (Bagla & 
Ray, 2005). Here D is the number of spatial dimensions in which 
gravitational clustering is being studied. If we consider an Einstein- 
deSitter universe then such a model is characterised by the scale 
of non-linearity tnl and the index n, where rjvz, is defined as the 
scale at which rms fluctuations in linearly extrapolated density con- 
trast (cri,(r)) is unity. 



TNL 



2 -(n + D) 

a r y ' 

2/(n+D) 



(1) 



If the evolution of clustering is self-similar then the non-linear cor- 
relation function plotted as a function of t/tnl should have the 
same shape at all epochs for a given index n. Note that the evo- 
lution is expected to be self-similar if tnl is the only scale in the 
problem and it is not expected in cosmologies with non-zero spatial 
curvature or a non-power-law expansion factor. 

Several studies have attempted to understand the nature of the 
asymptotic regime, mainly with help of N-Body simulations (Pea- 
cock & Dodds, 1994; Jain, Mo, & White, 1995; Colombi, Bouchet, 
& Hernquist, 1996; Padmanabhan et al., 1996; Peacock & Dodds, 
1996; Jain, 1997; Kanekar, 2000; Jing, 2001; Smith et al., 2003). 
Key conclusions of these studies about the asymptotic regime may 
be summarised as follows. 

• Stable clustering is not reached in the range of non-linearities 
explored by N-Body simulations. However, the departure from sta- 
ble clustering for most models is small. 

• Gravitational clustering does not erase memory of initial con- 
ditions, i.e., the non-linear power spectrum is not independent of 
the linear power spectrum to the extent it could be probed. There is 
no universal asymptotic regime. 

One of the reasons for the inability to resolve the problem of 
asymptotic regime has been the limited dynamic range of N-Body 
simulations. The problem of dynamic range can be circumvented 
by simulating a two-dimensional system (Bouchet, Adam, & Pel- 
lat, 1985; Engineer, Srinivasan, & Padmanabhan, 1999; Padmanab- 
han & Kanekar, 2000) instead of a three-dimensional one, wherein 
a much higher dynamic range can be achieved with similar compu- 
tational resources (Bagla, Engineer, & Padmanabhan, 1998; Mun- 
shi et al., 1998). Generic features like non-linear scaling relations 
are likely to be independent of dimension and — in fact — the 
scaling relations in the quasi-linear regime were predicted (Pad- 
manabhan, 1996) well before these could be tested in simulations 
(Bagla, Engineer, & Padmanabhan, 1998). If we can understand 
the nature of the asymptotic regime in two dimensions, it will help 
us solve the problem in three dimensions even if we cannot map 
the solution directly to the full problem in three dimensions. Previ- 
ous studies of gravitational clustering in two dimensions concluded 
that there is no stable clustering (Bagla, Engineer, & Padmanab- 
han, 1998; Munshi et al., 1998). Both the studies mentioned above 



were limited to £ 40 and the dynamic range in the non-linear 
regime was limited compared to the present study (see eqn.(6) for 
the definition of £). In this work we revisit the same issues, but 
using the 2D TreePM code (Ray, 2004) for N-body simulations. 
The TreePM code has a better force resolution as compared to a 
PM code (Bouchet, Adam, & Pellat, 1985), therefore in our simu- 
lations we have a significantly larger dynamic range over which we 
can study the asymptotic regime. 

A related issue is that of density profiles of haloes surrounding 
galaxies and clusters. Non-linear scaling relations in the asymp- 
totic regime depend on the kind of dynamical equilibrium that is 
reached in massive haloes, if an equilibrium is ever reached. Under 
the assumption of self-similar collapse, the density profiles of mas- 
sive haloes can be related to the initial conditions, and also to the 
dynamics within these haloes (Fillmore & Goldreich, 1984; Sub- 
ramanian, 2000). If these models are true then the slope of density 
profile at small scales is related to the index of power spectrum as: 

D (n + D + 1) 



p(r) oc r v ; rj 



(2) 



2 + D + n 

If self-similar collapse of haloes is the correct description then the 
final density profile of a halo carries information about the ini- 
tial profile and hence about the power spectrum (Bardeen et al., 
1986). It has been claimed that gravitational clustering in hierar- 
chical models leads to a universal density profile (Navarro, Frenk, 
& White, 1996) irrespective of the initial conditions, we shall refer 
to this density profile as the NFW profile. 

1 



r (r) oc 



{r/r s )(l + r/r s y 



(3) 



Here r s is a scale radius. We will also use the term r2oo for this. 
This profile is characterised by a r~ 3 decline at large radii and a 
cuspy inner profile of the form p(r) oc 1/r. The claim of univer- 
sality has been tested in several studies (Navarro, Frenk, & White, 
1996; Cole & Lacey, 1996; Fukushige & Makino, 1997; Tormen, 
Bouchet, & White, 1997; Moore et al., 1998; Kravtsov et al, 1998; 
Ghigna et al., 2000; Subramanian, Cen, & Ostriker, 2000; Klypin et 
al., 2001; Power et al., 2003; Fukushige, Kawai, & Makino, 2004; 
Navarro et al., 2004) and though there are disagreements, it appears 
that the NFW profile is consistent with the density profiles of mas- 
sive haloes in N-Body simulations even if it is not the best fit. This 
essentially implies that there is a scatter in density profiles and one 
can fit different functional forms to the N-Body data (Jing, 2000; 
Bullock et al., 2001). 

If true, the existence of NFW like profile implies that there are 
some universal aspects of gravitational clustering in an expanding 
universe. There is no ab initio derivation of the NFW profile, but 
one can argue that if there is a universal profile it should have same 
asymptotes as the form shown above (Syer & White, 1998; Nusser 
& Sheth, 1999; Padmanabhan, 2002a). 

We briefly review non-linear scaling relations in the following 
section. Numerical simulations used in the present work are dis- 
cussed in 83. We conclude in §4. 



2 NON-LINEAR SCALING RELATIONS 

In this section we review the non-linear scaling relations for grav- 
itational clustering in two dimensions. The non-linear and the lin- 
ear correlation functions at two different scales can be related by 
non-linear scaling relations (Hamilton et al., 1991; Nityananda & 
Padmanabhan, 1994). The relation between these scales is given by 
the characteristics of the pair conservation equation (Peebles, 1980; 
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Nityananda & Padmanabhan, 1994) under the assumption that the 
pair velocity function h (see below) depends on r and t through £ 
(Nityananda & Padmanabhan, 1994). For the two dimensional sys- 
tem of interest, this equation can be written as 



dD 
OA ' 

Here 



h(A,x)^=2h(A,X). 



£> = ]og(l + £), X = \og{x), A = log(a), 



(4) 



(5) 



and h — —v p /ax is the ratio between mean pair velocity and Hub- 
ble velocity. The volume averaged correlation function is defined 
as: 



= 2x~ 



/X 



)dr 



(6) 



where £; is the correlation function. The characteristics of this equa- 
tion are x 2 (l + ^(x,a)) = I 2 , where x and I are the two scales 
used in non-linear scaling relations. The self-similar models (Fill- 
more & Goldreich, 1984) imply that for collapse of cylindrical 
perturbations the turn-around radius and the initial density con- 
trast inside that shell are related as x ta oc I /Si oc Here 
£l is the linearly extrapolated mean correlation function. Noting 
that in two dimensions mass enclosed in a shell M oc i 2 , we find 
£(x) oc [fz,(0] m tne regime dominated by infall. Thus in 2-D 
the scaling relations are 



|i(o,I) 
h{a,l) 2 



(Linear) 
(Radial Infall) 
(Asymptotic limit) 



(7) 



The usual stable clustering limit is ft = 1 where the infall and 
Hubble expansion balance each other in structures in dynamical 
equilibrium. In such a case, or for any value of ft that does not 
depend on the initial conditions, slope of the non-linear correlation 
function is a unique function of the slope of the initial correlation 
function. For a general h we can relate the slope of the correlation 
function in the asymptotic regime to the slope of the initial linear 
correlation function 



£(a, a;) oc a; 



-2h(n+2) I (2 + h(n+2)) 



(8) 



If, however, h(n + 2) = constant then the non-linear correla- 
tion function has the same slope independent of the initial corre- 
lation function. This will happen if gravitational clustering erases 
all memory of the initial conditions. Note that this will make ft, and 
hence the velocity fields a function of the initial conditions. 

In our earlier work we had concluded that ft ~ 0.75 in the 
asymptotic regime (Bagla, Engineer, & Padmanabhan, 1998). In 
the same work we confirmed the prediction for the infall dominated 
quasi-linear regime, we found that the index in this regime is indeed 
2. (Note that in 3D, the indices for three regimes are 1, 3 and 3/i/2 
respectively.) In this work, we shall, among other things address 
the following two questions: 

• Does the asymptotic value of ft scales with n such that h(n + 
2) = constant? If it does, then the final slope of the non-linear 
correlation function will be independent of the initial slope. We 
will show that this does not happen. 

• Does ft reach a universal value independent of n? This seems 
to be the case and we will show that the behaviour of h is indepen- 
dent of n. But since all simulations, including this one probes only 
up to finite £, its actual asymptotic value is difficult to determine. 



3 NUMERICAL SIMULATIONS 

Our aim in this paper is to study the nature of scaling relations in 
the asymptotic regime. We choose to work in two dimensions in 
order to improve the dynamic range over which the results of N- 
Body simulations are reliable. We wish to do this without changing 
the system completely as the ultimate goal is to carry over the un- 
derstanding to three dimensional systems. Therefore we study the 
evolution of two dimensional perturbations in a three dimensional 
expanding universe (Fillmore & Goldreich, 1984). We start with a 
set of infinitely long straight "needles" all pointing along one axis. 
The evolution keeps the "needles" pointed in the same direction and 
we study gravitational clustering of these needles in an orthogonal 
plane. Particles in the N-body simulation represent the intersection 
of these "needles" with this plane and the force of interaction be- 
tween "particles" is given by the solution of the Poisson equation 
in two dimensions. The inter-particle force falls as 1 jr, unlike the 
inverse square force in three dimensions. 

We use the two dimensional TreePM code (Ray, 2004), a mod- 
ified version of the three dimensional TreePM code (Bagla, 2002; 
Bagla & Ray, 2003). The TreePM code splits the interaction force 
into two parts, the long range force is computed in Fourier space 
and the short range component is computed in real space. The 
TreePM code does not use the usual PM force as the long range 
force, this allows us to control errors in force over the entire range 
of scales (Bagla, 2002; Bagla & Ray, 2003; Ray, 2004). 

We soften the two-dimensional gravitational force at small 
scales in order to ensure collisionless evolution of the particle distri- 
bution in our simulations. We use a cubic spline form for softening, 
i.e., each particle is assumed to have an extended mass distribution 
represented by the normalised spline kernel used in the SPH for- 
malism (Monaghan, 1992). This softens the force at scales smaller 
than the softening length e, and at scales larger than e we get the 
full 1 jr force. All N-Body simulations reported in this paper used 
e > 0.2 in units of the average inter-particle separation. 

We have carried out a series of numerical simulations with 
power law initial spectra with indices n = —0.4, 0.0 and 1.0. Sev- 
eral realisations of each power spectrum were simulated. Scatter in 
values across this suite of simulations was used to estimate the error 
bars. It is difficult to reach the asymptotic regime forn <C —0.4 be- 
fore the perturbations at the box scale become important and hence 
we have not used models with a more negative index. Models were 
normalised so that A 2 (fc = 8kf,a = 1) = 1 where kf is the 
wave number of the fundamental mode, essentially the fluctuations 
at 1/8 of the box size were normalised to unity at a — 1. Simula- 
tions were run upto a sufficiently late epoch while requiring that the 
fluctuations at half the box scale were well within the linear regime 
at the final epoch, i.e., A (k — 2k f, a,f in ) <C 1. 

Simulations discussed here used 2048 2 particles on a 2048 2 
grid in an Einstein de Sitter background universe. TreePM code 
was used for all the simulations, though we have checked many of 
the results in larger 4096 PM simulations. We used r a — 1.0 grid 
length in the TreePM code, this is the scale used for dividing the 
force into a short range and a long range component (Bagla, 2002; 
Bagla & Ray, 2003; Ray, 2004). Error in force depends largely on 
the choice of this scale and the cell opening criterion, we use 6 C — 
0.5. 

In an Einstein-deSitter universe the evolution of perturbations 
with a power law power spectrum is expected to be scale-invariant. 
The only scale introduced by gravitational clustering is the scale of 
non-linearity rjvi. 
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r / r NL 



Figure 1. f is plotted as a function of r/r NL for n = —0.4 (upper panel) 
and n = 1.0 (lower panel) simulations. These curves span more than a 
decade in scale factor. Clearly, over the entire range of clustering the evolu- 
tion is self-similar. Dashed lines in these panels shows the asymptotic slope 
for h = 0.75, a result suggested by an earlier study. 



We outline the key results of these simulation studies in the 
following subsections. 



3.1 Correlation Function 

We find that the evolution of the system is self-similar in that 
all the relevant quantities have the same form when scaled by 
r NL oc a 2l( - n+2 \ where r NL is the scale where the variance of 
density fluctuations is unity. Fig. 1 shows the averaged correla- 
tion function £(r) as a function of r/r NL for simulations with 
n = —0.4 (upper panel) and n = 1.0 (lower panel). Curves were 
plotted for scales larger than twice the softening length in order to 



keep out features that depend on the choice of softening length. The 
evolution is self-similar up to and beyond £ = 100. 

The shape of the correlation function f(r) for weaker non- 
linearities is consistent with the results of earlier studies (Bagla, 
Engineer, & Padmanabhan, 1998; Munshi et al., 1998). There was 
disagreement between the two studies for the asymptotic regime: 
one study (Bagla, Engineer, & Padmanabhan, 1998) favoured h ~ 
0.75 and the other found h(n + 2) = 1.3 (Munshi et al., 1998). The 
latter behaviour erases memory of initial conditions and the slope 
of the correlation function in the asymptotic regime is the same for 
all models. The dynamic range in both the studies was limited to 
f~10. 

The current work with higher dynamic range agrees with the 
slope for h = 0.75 in the region of overlap with the previous work 
(Bagla, Engineer, & Padmanabhan, 1998) and rules out the pos- 
sibility of h(n + 2) = 1.3. The dashed line with the expected 
slope for h = 0.75 is marked in Fig. 1 and it runs parallel to the 
curve up to about £ ~ 40, the largest non-linearity studied ear- 
lier (Bagla, Engineer, & Padmanabhan, 1998). As clustering in- 
creases, the slope of the correlation function decreases below the 
slope expected for h = 0.75. Slope of the correlation function for 
these two models is shown in Fig. 2. The asymptotic slope of the 
correlation function is —0.53 ^ 7 ^ —0.50 for n = —0.4 and 
—0.80 ^ 7 ^ —0.77 for n = 1, where 7 = 9 log £/d log r and is 
evaluated at £ ^ 100. We have given the 95% confidence limits for 
the slope, the limits were derived using a \ 2 At to the correlation 
function in the asymptotic regime. Different values of 7 imply de- 
parture from the h(n + 2) = constant asymptote. The values of 7 
listed above are consistent with 0.416 < h < 0.451 for n — —0.4 
and 0.417 < h < 0.444 for n = 1. Thus we recover a similar 
range of asymptotic values for h from the correlation function. As 
we will show in §3.2, the slope of the correlation function is a more 
stable estimator of h than a direct determination. 

Fig. 2 shows the non-linear scaling relations. £(x, a) is plot- 
ted as a function of a) for n = —0.4, n = and n = 1. 
Data from multiple epochs have been used here. Curves marking 
the linear, quasi-linear and asymptotic regime are shown here. We 
have also marked a line showing the stable clustering limit. The 
equations for the piecewise fit are: 

C a(o,i); £(a:)<0.5 
|(o, as) = I 2 ^(o, if; 0.5 «S f(a;) < 8 (9) 

[ 10.5 & (a,/) 045 ; 22sC<e» 

It is remarkable that the points for all the three models follow the 
same non-linear scaling relation. The asymptotic slope implied by 
the scaling relation lies near h = 0.45 but it is difficult to say 
whether this is the final value or if this continues to decrease as 
we go to higher non-linearities. This evolution is inconsistent with 
the asymptote of h(n + 2) = constant, as the scaling relations for 
different models would have been different in that case. 

The results shown here with the TreePM code match with the 
correlation function obtained in 4096 2 PM simulations. This com- 
parison has been carried out up to £ < 70. 



3.2 Pair Velocity 

Fig. 3 shows the pair velocity h(r, a) as a function of f(r, a) for a 
large number of epochs a for n = —0.4 (upper panel) and n = 1.0 
(lower panel). The dashed line shows the expected value of h in 
the linear limit. For the linear and quasi-linear regime (f < 10) 
ft is a single valued function of f, irrespective of the epoch for a 
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Figure 2. The upper panel shows the non-linear scaling relation. £(x, a) 
is plotted as a function of §l('> a) for n = —0.4, n = and n = 1. 
Data from multiple epochs has been used here. Curves marking the lin- 
ear (thick curve), quasi-linear (dashed curve) and asymptotic regime (dot- 
dashed curve) are shown here. We have also marked a dotted line showing 
the stable clustering limit. A remarkable feature of clustering in two dimen- 
sions is that the non-linear correlation function in the asymptotic regime 
is smaller than the linear correlation function. The slope 7 is shown as a 
function of £(x, a) in the lower panel. Circles mark the slope for the n = 1 
model and the plus sign marks the slope for the n = —0.4 model. Points 
for the n = 1 model are displaced upwards by 2 for clarity. 



given model. Indeed, we find that there is no significant difference 
between the curves for different models at £ < 10, and at higher 
non-linearities there is a large scatter in the curves so it is diffi- 
cult to test any claims using Fig. 3. Thus simulations are consis- 
tent with the ansatz that h depends on epoch and scale through £ 
(Nityananda & Padmanabhan, 1994) and this allows us to find the 
form of the non-linear scaling relation. We can certainly do this for 




0.1 1 10 100 

f(r) 




F(r) 



Figure 3. The pair velocity h(r, a) is plotted as a function of f(r, a) for a 
large number of epochs a for n = —0.4 (upper panel) and n = 1.0 (lower 
panel). The dashed line shows the expected value of h in the linear limit. 
The h — £ curves are the same across all epochs and for all models within 
the scatter. 



a given model and as we do not find significant differences between 
different models, there does not seem to be any power spectrum de- 
pendence in scaling relations. However, there will be a dependence 
on cosmology (Padmanabhan et al., 1996). Conversely, non-linear 
scaling relations can be used to find h(£). 

We have checked that the scaling relations plotted in Fig. 2 are 
consistent with the curves plotted in Fig. 3. This is demon- 
strated by Fig. 4 where h determined directly is compared with 
values obtained from the non-linear scaling relation. 

The value of h fluctuates a lot at late epochs, we mentioned in 
the previous section that the slope of the correlation function gives 
a more stable estimate of h. Due to this, some of the late epochs 
used in Fig. 1 have not been plotted in Fig. 3. The key feature in 
this regime is that h < 1 and it does not show any sign of heading 
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100 



Figure 4. The pair velocity h(r, a) is plotted as a function of a) at late 
times for n = —0.4. Error bars are also shown here. The range of h(r, a) 
expected from eqn.(4) is shown here using dashed lines. The expected val- 
ues and the directly measured values of h are consistent within the error 
bars for nearly the entire range of scales. 



towards the stable clustering limit of h = 1. This is consistent with 
the results from the slope of the correlation function. 

It is important to understand what h < 1 means in terms of 
dynamics. One possibility is that the haloes are evaporating. This 
can happen if two-body relaxation is important, but we have used 
e ^ 0.2 grid lengths in all our simulations and number density of 
particles in most haloes is very high. Hence, two-body scattering 
should not be important. A poor integrator for the equation of mo- 
tion can also lead to the evaporation of haloes. We have tested the 
integrator in three-body problems and with highly eccentric binary 
orbits to ensure that the problem does not lie with incorrect evolu- 
tion of trajectories. 

In order to probe the dynamics further, we also study non- 
radial motions. Analogous to h, the radial component of the pair 
velocity scaled by Hubble velocity, we define the transverse pair 
velocity function as: 



g(r,a) = 



cVp V Vy X v) 



Hr 2 



(10) 



Here the angle brackets denote averaging over all pairs with sep- 
aration r, v p is the pair velocity (peculiar velocity) and r p is the 
separation of the pair of particles. H is the Hubble parameter. The 
equivalent object in three dimensions is: 



g(r,a) = |g| = 



Hr 2 



(11) 



In two dimensions the "cross" product of two vectors is just a num- 
ber. Fig. 5 shows \g\ as a function of £. We have plotted the magni- 
tude of g(r, a) as its value oscillates about zero at large-scales. This 
curve is plotted for only one epoch but the relative value of \g\ com- 
pared to h shows that this is a significant component. Note that we 
require coherent transverse motions in order to detect anything here 
as random transverse motions will cancel out in the sum. We do 
not find any systematic excess in pairwise transverse velocity dis- 
persion as compared to the pairwise collinear velocity dispersion. 




Figure 5. The transverse component of pair velocity g(r, a) is plotted as 
a function of §(r, a) for a late epoch for n = —0.4 (upper panel) and 
n = 1.0 (lower panel) models. The magnitude \g\ is plotted here, a remark- 
able fact apparent from this figure is that the transverse component is more 
important for n = 1. At £ 2> 1, the magnitude of g is comparable with that 
ofh. 



In dynamical equilibrium these should have the same magnitude in 
two dimensions and at £ ^ 10 we find this to be true. 

To investigate this further, we have studied motions in individ- 
ual haloes by plotting the angular momentum profile as a function 
of distance from the centre of the halo. There are no systematic fea- 
tures of shape worth commenting except that for most haloes there 
are annuli with coherent rotation. Note that this does not imply a 
coherent rotation for the halo as the direction of rotation reverses a 
few times. Fig. 6 shows circularly averaged angular momentum of 
particles in a halo within a distance r from the centre of mass of the 
halo for some (randomly selected) haloes. The upper panel is for in- 
dex n = —0.4 and the lower for n = 1.0. Clearly, all the haloes 
in both panels have annuli with coherent rotation. In some cases, 
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Figure 6. In this figure, we plot circularly averaged angular momentum of 
particles in a halo at a distance r from the centre of mass of the halo for 
some haloes. Distance r is scaled by r200 an d the contribution of particles 
at a distance r to the angular momentum is scaled by the total angular mo- 
mentum contained within r200- The upper panel is for index n = —0.4 
and the lower for n = 1.0. Different line styles here correspond to different 
haloes. 

angular momentum contributed by adjacent annuli is in opposite 
direction. The net angular momentum in these haloes contributes a 
small fraction of the kinetic energy of particles in these haloes. A 
discussion of the selection criteria used to identify haloes is given 
in the next sub-section. 

It is important to understand the origin of this apparent coher- 
ent rotation in massive haloes. Comparison of g — £ plot (Fig. 5) 
for n = —0.4 and n = 1 offers us a clue, transverse motions are 
much stronger for n = 1 as compared to the other model. Col- 
lisions between sub-structure falling into haloes in the early phase 
can generate such transverse motions (Bagla, Prasad, & Ray, 2005), 
and there is certainly more substructure in the n = 1 model. We 
put forward the hypothesis that collisions between substructure are 
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Figure 7. Spherically averaged density profiles of a number of spherically 
symmetric virialised haloes from simulations with indices n = —0.4 (up- 
per panel) and n = 1.0 (lower panel). Density is plotted as a function of 
r /V200 • The four dashed lines in each panel mark the approximate extremes 
of inner (rf) and outer slopes (/3) for the density profiles. The smallest haloes 
have more than 10^ particles inside the virial radius, and the largest haloes 
have more than 10 4 particles. All the haloes in the figure were taken from 
TreePM simulations with 2048 2 particles. 

responsible for generation of coherent transverse motions. We are 
testing this hypothesis in a series of numerical simulations. Results 
of these studies will be reported elsewhere. 

3.3 Density Profiles of Massive Haloes 

Density profiles of massive haloes and non-linear scaling relations 
are closely related. Density profiles of very massive haloes are eas- 
ier to study as these extend over a large region, large compared to 
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the softening length which is the smallest scale we can resolve in 
our simulations. We identified haloes in our simulations using the 
method of density peaks (Tormen, Bouchet, & White, 1997). In this 
method the density field is smoothed at a large length scale R = 5 
grid lengths and peaks in the smoothed density field are identified. 
Particles within a distance R of a peak are selected and we find 
the centre of mass for these particles. We now select a smaller set 
of particles out of these, those that are within R — AR from the 
centre of mass and compute the centre of mass. These iterations 
continue till we are left with only a small number of particles, say 
about 100. The final centre of mass is taken to be the centre of the 
halo and the density profiles are plotted as a function of distance 
from this centre. We apply two further criteria: 

• Central density contrast of the halo must be very large, S > 
500. 

• Haloes and their neighbourhood should be smooth and there 
should be no major merger going on. To implement this in an ob- 
jective manner we compute the ratio \ = (cos 2 <j>)/(sm 2 <j>) where 
average is over all particles within the central region of the halo and 
4> is the position angle of a particle from the centre of mass. We re- 
quire 0.9 ^ x ^ 1-1 f° r a halo t0 be used for density profiles. 

These two criteria reject about 75% of all the haloes with rejection 
being more frequent for models with a larger index n. 

Fig. 7 shows the density profiles of haloes that satisfy these 
criteria in our simulations. Density profiles of a large number of 
haloes are plotted, for each halo we have drawn density p as a func- 
tion of r/r2oo, where V2oo is the scale at which the density is 200 
times the average. Average density is unity in these simulations. It 
is clear that there is a large scatter in density profiles. The inner 
slope varies from 0.4 to 1.4 for n = —0.4 whereas it is in the 
range 0.6 to 1.6 for n = 1. The inner slope ij is obtained by fit- 
ting a power law that passes through p = 200 and r/rwo = 1-0 
and only points inside of this radius are used. The distribution of 
the inner slope of density profiles is very broad, as is obvious from 
the range. The median r\ is approximately 0.9 for n — —0.4 and 
1.1 for n — 1.0. Self-similar, spherical collapse predicts an inner 
slope of r] — 0.89 and r\ — 1.2 for n = —0.4 and n = 1, respec- 
tively (eqn.(2)). These values are well within the scatter of inner 
slopes. Clearly, the answer lies in between the prediction for the 
self-similar collapse and a universal slope. 

It is difficult to explain the large scatter seen in density pro- 
files. The origin of scatter can lie in environmental dependence as 
well as in different formation histories for haloes. We can, however, 
address the issue of the small systematic shift in median slope as 
compared to that expected from stable clustering. 

Transverse motions can affect the density profile of massive 
haloes in a significant manner (Subramanian, 2000). Such effects 
will be stronger if the transverse motions are stronger, and as is 
obvious from the Fig. 5 transverse motions are stronger for larger 
n. Rotation is playing a more important role in models with more 
substructure and is leading to flatter density profiles. This reduces 
the difference between density profiles of massive haloes for these 
models. This, we believe, is a common explanation for the non- 
linear scaling relations and density profiles. 



3.4 Virial Equilibrium 

The stable clustering limit requires formation of objects that are in 
dynamical equilibrium. Following Peebles (1980), we compute the 
second time derivative of the moment of inertia and set it to zero 




Figure 8. The average kinetic energy T of massive haloes is scaled by 
G M 2 and plotted against M 2 on a log-log scale. Each point on this graph 
represents one halo, we have used the set of haloes used for studying den- 
sity profiles. This plot shows that there is large scatter about the average 
relation T <x M 2 . 



for a cluster of particles in dynamical equilibrium. This gives us the 
condition to be satisfied for dynamical equilibrium. We get: 



(12) 



where the sums are over all particles and dots denote time deriva- 
tives. Particles are located at and have velocities and rrii is the 
mass of the ith particle. Thus the kinetic energy depends only on 
the total mass of the halo. If the massive haloes that we are study- 
ing here are in virial equilibrium then the kinetic energy T for these 
haloes should be proportional to M 2 . For each halo used for density 
profiles, all particles within the radius 7-200 were used to plot Fig. 8. 
It is clear that these quantities are correlated even though there is a 
large scatter about the average relation. Thus we may conclude that 
the massive haloes are close to dynamical equilibrium. 



4 CONCLUSIONS 

The basic motivation for this work was to improve our understand- 
ing of non-linear scaling relations for two dimensional gravita- 
tional clustering. Previous studies (Bagla, Engineer, & Padmanab- 
han, 1998; Munshi et al., 1998) established the behaviour for the 
quasi-linear regime, but could not probe the nature of the asymp- 
totic regime. We used the 2D TreePM code (Ray, 2004) which al- 
lows us to study the non-linear regime in far greater detail, as com- 
pared to the PM code used in earlier studies. The results of our 
study can be summarised as follows: 

• The evolution of the correlation function for power-law mod- 
els is scale invariant up to the highest non-linearities reached in our 
simulations. 

• We reproduce the results of earlier studies in the quasi-linear 
regime. We confirm that the slope of the non-linear scaling relation 
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in this regime is close to 2, as predicted on the basis of the infall 
dominated growth model (Padmanabhan, 1996). 

• We do not find any difference in the non-linear scaling relation 
for different power law models. Thus the information about initial 
conditions is retained even in the extremely non-linear regime. 

• The stable clustering limit (h — * 1 as £ — * oo) is not reached. 
(Though it is generally claimed that stable clustering with h — ► 1 is 
a natural asymptotic state, it has been shown that the Davis-Peebles 
scale-invariant solution (Davis & Peebles, 1977) and the hierarchi- 
cal model for the three-point function are inconsistent with the stan- 
dard stable-clustering picture (Kanekar, 2000). This, of course, is 
for gravitational collapse in three dimensions.) This can be con- 
cluded from the slope of the correlation function, pair velocity, as 
well as density profiles of massive haloes. 

• Pair velocities in the asymptotic regime are smaller than ex- 
pected from the stable clustering model. This does not imply that 
massive haloes are evaporating. 

• The transverse component of pair velocity is more significant 
for models with a larger spectral index, i.e., for models with more 
substructure. This suggests that gravitational collisions between 
substructure orbiting within haloes are responsible for generating 
coherent transverse motions. 

• We find that there is no universal density profile for massive 
haloes in two dimensional gravitational clustering. There is a large 
scatter in inner as well as outer slopes of density profiles. Median 
value of the inner slope is different for different models, though the 
difference is much smaller than the scatter in values. Tests show 
that these clusters are close to dynamical equilibrium. 

• The difference in the median inner slope for different models 
is smaller than expected from self-similar spherical collapse. We 
argue that the difference in significance of transverse motions for 
different models is responsible for this change. 

Many of the conclusions listed here have been checked inde- 
pendently in large PM simulations (4096 2 particles) up to £ ^ 100. 

Further work is required to test our hypothesis that gravita- 
tional collisions between substructure are responsible for gener- 
ating coherent transverse motions. We are carrying out a detailed 
study using numerical simulations where we tune the amount of 
substructure to test this proposal. 

Considerable work has been done to study the non-linear 
scaling relations in three dimensions (Hamilton et al., 1991; 
Nityananda & Padmanabhan, 1994; Peacock & Dodds, 1994; Jain, 
Mo, & White, 1995; Peacock & Dodds, 1996; Padmanabhan, 1996; 
Padmanabhan et al., 1996; Bagla, Engineer, & Padmanabhan, 1998; 
Smith et al., 2003). We intend to study these again in light of insight 
developed using tests in clustering in two dimensions. It has been 
shown that the scaling relations for gravitational clustering in three 
dimensions depend on the initial spectrum of fluctuations (Smith 
et al., 2003). It is important to understand the origin of this index 
dependence and test whether substructure plays a role in this. The 
transverse pair velocity function has been introduced in eqn.(10), 
we are studying the three dimensional equivalent of this quantity in 
order to probe the role of transverse motions in gravitational clus- 
tering in the full three dimensional problem. These studies should 
help us understand the non-linear scaling relations better. 
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